% Example83.m
% Implements JPEG2000 wavelet transform-based compression.
% Implements only the scalar quantization normative of the
% standard. Implements up to NL levels of 2D DWT using
% Daubechies 9/7 filter bank.
% Implements 4:2:0 and 4:2:2 sampling formats specified
% by the user.
% Uses the frequency weightings recommended by JPEG2000 standard
% in the uniform scalar quantization of the DWT coefficients.
%
clear
NL = 3; % number of levels of 2D DWT
SamplingFormat = '4:2:2';
A = imread('cameraman.tif');
%A = imread('lighthouse.ras');
% make sure that the image size is divisible by NL
[X,Y,Z] = size(A);
if mod(X , 2^NL)~=0
Height = floor(X/(2^NL))*(2^NL);
else
Height = X;
end
if mod(Y,2^NL)~=0
Width = floor(Y/(2^NL))*(2^NL);
else
Width = Y;
end
Depth = Z;
clear X Y Z
% Do forward Irreversible Component Transformation (ICT)
% In JPEG2000 jargon, y0 corresponds to the Luma Y,
% y1 to Cb and y2 to Cr.
if Depth == 1
y0 = double(A(1:Height,1:Width)) - 128; % DC level shift
else
A1 = double(A(1:Height,1:Width,:))-128; % DC level shift
y0 = 0.299*A1(:,:,1) + 0.587*A1(:,:,2) + 0.144*A1(:,:,3);
y1 = -0.16875*A1(:,:,1) - 0.33126*A1(:,:,2) + 0.5*A1(:,:,3);
y2 = 0.5*A1(:,:,1) - 0.41869*A1(:,:,2) - 0.08131*A1(:,:,3);
clear A1
switch SamplingFormat
case '4:2:0'
y1 = imresize(y1,[Height/2 Width/2],'cubic');
y2 = imresize(y2,[Height/2 Width/2],'cubic');
case '4:2:2'
y1 = imresize(y1,[Height Width/2],'cubic');
y2 = imresize(y2,[Height Width/2],'cubic');
end
end
% Daubechies 9/7 Filters
% Note: The Daubechies 9/7 filter bank has 9 taps for its analysis
% Lowpass and synthesis highpass FIR filters;
% 7 taps for its analysis highpass and synthesis lowpass
% FIR filters. However, in order to obtain subbands that are
% exact multiples of 2^-n, n being a positive integer, the
% number of taps in the filters must be even. That is why you
% will notice the filter lengths to be 10 instead of 9 or 7.
% The extra tap values are zero.
%
%LOD = analysis lowpass filter
%LOR = synthesis lowpass filter
%HID = analysis highpass filter
%HIR = synthesis highpass filter
LOD = [0 0.0267487574 -0.0168641184 -0.0782232665 0.2668641184...
0.6029490182 0.2668641184 -0.0782232665 -0.0168641184...
0.0267487574];
HID = [0 0.0912717631 -0.0575435263 -0.5912717631...
1.1150870525 -0.5912717631 -0.0575435263 0.0912717631 0 0];
%
LOR = [0 -0.0912717631 -0.0575435263 0.5912717631...
1.1150870525 0.5912717631 -0.0575435263 -0.0912717631 0 0];
HIR = [0 0.0267487574 0.0168641184 -0.0782232665 -0.2668641184...
0.6029490182 -0.2668641184 -0.0782232665 0.0168641184...
0.0267487574];
%
% "dwtCoef" cell array has NL x 4 cells.
% In each cell at level NL, the 1st component is the LL
% coefficients,
% 2nd component is HL coefficients, 3rd component LH and
% 4th component HH coefficients.
%
% Do an NL-level 2D DWT of the Y component
dwty0 = cell(NL,4);
for n = 1:NL
if n==1
[dwty0{n,1},dwty0{n,2},dwty0{n,3},dwty0{n,4}] =...
dwt2(y0,LOD,HID,'mode','per');
else
[dwty0{n,1},dwty0{n,2},dwty0{n,3},dwty0{n,4}] =...
dwt2(dwty0{n-1,1},LOD,HID,'mode','per');
end
end
% Quantize the Y DWT coefficients
[dwty0,y0Bits] = jpg2KQuantize(dwty0,NL,1);
%
% Do an NL-level 2D IDWT of the Y component
% The cell array idwtimg has NL cells with each cell
% corresponding to
% the reconstructed LL image at that level.
idwty0 = cell(NL,1);
for n = NL:-1:1
if n==3
idwty0{n}= idwt2(dwty0{n,1},dwty0{n,2},dwty0{n,3},...
dwty0{n,4},LOR,HIR,'mode','per');
else
idwty0{n}= idwt2(idwty0{n+1},dwty0{n,2},dwty0{n,3},...
dwty0{n,4},LOR,HIR,'mode','per');
end
end
SNRy0 = 20*log10(std2(y0)/std2(y0-idwty0{1}));
if Depth > 1
% Do an NL-level 2D DWT of the Cb and Cr components
dwty1 = cell(NL,4);
dwty2 = cell(NL,4);
for n = 1:NL
if n==1
[dwty1{n,1},dwty1{n,2},dwty1{n,3},dwty1{n,4}] =...
dwt2(y1,LOD,HID,'mode','per');
[dwty2{n,1},dwty2{n,2},dwty2{n,3},dwty2{n,4}] =...
dwt2(y2,LOD,HID,'mode','per');
else
[dwty1{n,1},dwty1{n,2},dwty1{n,3},dwty1{n,4}] =...
dwt2(dwty1{n-1,1},LOD,HID,'mode','per');
[dwty2{n,1},dwty2{n,2},dwty2{n,3},dwty2{n,4}] =...
dwt2(dwty2{n-1,1},LOD,HID,'mode','per');
end
end
%
% Quantize the Cb and Cr DWT coefficients
[dwty1,y1Bits] = jpg2KQuantize(dwty1,NL,2);
[dwty2,y2Bits] = jpg2KQuantize(dwty2,NL,3);
%
% % Do an NL-level 2D IDWT of the Cb and Cr components
idwty1 = cell(NL,1);
idwty2 = cell(NL,1);
%
for n = NL:-1:1
if n==3
idwty1{n}= idwt2(dwty1{n,1},dwty1{n,2},dwty1{n,3},...
dwty1{n,4},LOR,HIR,'mode','per');
idwty2{n}= idwt2(dwty2{n,1},dwty2{n,2},dwty2{n,3},...
dwty2{n,4},LOR,HIR,'mode','per');
else
idwty1{n}= idwt2(idwty1{n+1},dwty1{n,2},dwty1{n,3},...
dwty1{n,4},LOR,HIR,'mode','per');
idwty2{n}= idwt2(idwty2{n+1},dwty2{n,2},dwty2{n,3},...
dwty2{n,4},LOR,HIR,'mode','per');
end
end
% Compute the SNR of Cb and Cr components due to quantization
SNRy1 = 20*log10(std2(y1)/std2(y1-idwty1{1}));
SNRy2 = 20*log10(std2(y2)/std2(y2-idwty2{1}));
% Upsample Cb & Cr components to original full size
y1 = imresize(idwty1{1},[Height Width],'cubic');
y2 = imresize(idwty2{1},[Height Width],'cubic');
end
% Do inverse ICT & level shift
if Depth == 1
Ahat = idwty0{1}+ 128;

subplot(1,2,1), imshow(A), title('Originnal Image');
subplot(1,2,2), imshow(uint8(Ahat)), title('JPEG2000 Compressed Image');

sprintf('Average bit rate = %5.2f bpp\n',y0Bits/(Height*Width))
sprintf('SNR(Y) = %4.2f dB\n',SNRy0)
else
Ahat = zeros(Height,Width,Depth);
Ahat(:,:,1) = idwty0{1}+ 1.402*y2 + 128;
Ahat(:,:,2) = idwty0{1}- 0.34413*y1 - 0.71414*y2 + 128;
Ahat(:,:,3) = idwty0{1}+ 1.772*y1 + 128;

subplot(1,2,1), imshow(A), title('Originnal Image');
subplot(1,2,2), imshow(uint8(Ahat)), title('JPEG2000 Compressed Image');

sprintf('Average bit rate = %5.2f bpp\n',...
(y0Bits+y1Bits+y2Bits)/(Height*Width))
SNRy0 = 20*log10(std2(y0)/std2(y0-idwty0{1}));
sprintf('SNR(Y) %4.2fdB\tSNR(Cb) = %4.2f dB\tSNR(Cr) =%4.2fdB\n',SNRy0,SNRy1,SNRy2)
end